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Abstract 


We  present  a  method  to  reconstruct  images  from  finite  sets  of  noisy  projections  that  may  be  avail¬ 
able  only  over  limited  or  sparse  angles.  The  algorithm  calculates  the  maximum  a  posteriori  (MAP) 
estimate  of  the  full  sinogram  (which  is  an  image  of  the  2-D  Radon  transform  of  the  object)  from 
the  available  data.  It  is  implemented  using  a  primal-dual  constrained  optimization  procedure  that 
solves  a  partial  differential  equation  in  the  primal  phase  with  an  efficient  local  relaxation  algorithm 
and  uses  a  simple  Lagrange  multiplier  update  in  the  dual  phase.  The  sinogram  prior  probability  is 
given  by  a  Markov  random  field  (MRF)  that  includes  information  about  the  mass,  center  of  mass, 
and  convex  hull  of  the  object,  and  about  the  smoothness,  fundamental  constraints,  and  periodicity 
of  the  2-D  Radon  transform.  The  object  is  reconstructed  using  convolution  backprojection  applied 
to  the  estimated  sinogram.  We  show  severed  reconstructed  objects  which  are  obtained  from  simu¬ 
lated  limited-angle  and  sparse-angle  data  using  the  described  algorithm,  and  compare  these  results 
to  images  obtained  using  convolution  backprojection  directly. 
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Research  Office  grant  DAAL03-86-K-1071.  In  addition,  the  work  of  the  first  author  was  partially  supported 
by  a  U.S.  Army  Research  Office  Fellowship. 
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1  INTRODUCTION 


This  paper  eiddresses  the  problem  of  reconstruction  from  projections,  the  theoretical  and 
practical  aspects  of  which  have  received  much  attention  over  the  past  two  decades.  Among 
the  applications  that  use  the  currently  available  reconstruction  techniques  are  medicine, 
optics,  material  science,  astronomy,  geophysics,  and  magnetic  resonance  imaging  [1].  The 
most  widely  known  application  of  this  theory  is  the  problem  of  medical  transmission  X- 
ray  tomography  [2].  In  this  discipline,  “pencil  beam”  X-rays  are  fired  from  many  angles 
through  a  single  cross  section  of  the  body,  effectively  measuring  line  integrals  of  the  2-D 
X-ray  density  function  corresponding  to  the  various  tissues  in  the  cross  section.  A  collection 
of  line  integreJs  obt£iined  over  lines  with  the  same  angle  but  different  lateral  positions  forms 
a  1-D  function  called  a  projection.  Given  a  set  of  projections  taken  from  many  different 
angles,  an  image  of  the  density  function  may  be  reconstructed  and  used  in  diagnosis.  For 
many  medical  conditions  this  tomographic  approach  to  imaging  of  the  body  leads  to  greatly 
improved  imagery  over  conventional  (chest-type)  X-ray  images  and  has  proven  to  be  of  great 
benefit  in  medical  diagnosis  [3]. 

Consider  a  function  f(x)  defined  on  the  plane  as  depicted  in  Fig.  1.  We  denote  the 
integral  of  /  along  the  line  L(t,  0)  by  g{t,  6).  The  function  g  for  all  values  of  t  and  $  is  called 
the  2-D  Radon  transform  of  /,  and  an  image  of  g{t,9),  with  t  in  the  vertical  direction  and 
9  in  the  horizontal  direction,  is  called  a  sinogram.  For  a  single  value  of  d,  ^  is  a  function  of 
t  only,  and  is  called  the  projection  of  /  at  angle  9.  The  Radon  transform  of  f{x)  may  be 
written  as 

g{t,9)=  f  f{x)6{t- x-oj)dx  (1) 

JeeiR* 

where  oj  =  [cos^  sintf]'^,  and  5(*)  is  the  Dirac  delta  function.  It  turns  out  that  only  certain 
functions  g(t,9)  are  valid  Radon  transforms;  there  are  inherent  mathematical  consistency 
conditions  that  constrain  g{t,9)  to  lie  in  a  particular  functional  subspace  defined  on  the 
cylinder  IR^  x  S^. 

The  fact  that  (1)  is  invertible  (for  a  wide  class  of  functions)  is  well  known.  Deans 
[4]  describes  many  of  the  known  exact  inversion  formulas.  Except  under  certain  (usually 
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impractical)  circumstances,  however,  it  is  not  possible  to  determine  /  exactly  given  the 
value  of  g  for  only  a  finite  number  of  lines  It  has  been  the  primary  concern  of 

engineers  and  physicists  in  this  field,  over  the  last  20  years  or  so,  to  study  approximate 
inversion  algorithms  given  such  a  finite  measurement  set.  The  performance  of  any  particular 
algorithm  depends  on  the  nature  of  the  measurements  -  their  number  and  arrangement, 
and  their  noise  properties  -  and  often  on  the  nature  of  the  object  itself. 

In  this  paper,  we  are  concerned  with  the  case  of  low  signed  to  noise  ratio  (SNR)  and 
limited-angle  or  sparse-angle  measurement  configurations  with  parallel-ray  projections.  In 
the  medical  CT  problem,  for  example,  a  line  integral  measurement  may  be  noisy  if  low 
energy  X-rays  are  used.  Data  acquisition  may  be  restricted  to  a  limited  angular  range  if 
there  is  an  obstruction,  for  example,  and  may  be  sparsely  sampled  if  there  is  a  requirement 
to  obtain  the  data  extremely  rapidly.  In  these  situations,  images  reconstructed  using  con¬ 
ventional  methods  are  degraded  by  a  variety  of  artifacts  [5],  and  alternate  methods  must 
be  used. 

Several  investigators  have  developed  edgorithms  to  compensate  for  some  of  the  data 
deficiencies  described  above.  The  modified  transform  methods  of  [6],  [7],  emd  [8]  ziccount 
for  missing  projections  but  do  not  explicitly  address  the  the  issue  of  noise  in  the  data. 
Finite-series  expansion  methods  use  additional  criteria  such  as  minimum  norm  [9],  minimum 
variance  [10,11],  and  maximum  entropy  [12]  to  account  for  noise,  or  missing  projections,  or 
both.  In  many  cases,  streaking  artifacts  are  still  present  in  the  reconstructed  images  [13] 
and,  in  all  but  certain  very  special  imaging  geometries,  such  as  in  [14],  the  computations 
are  very  time  consuming  and,  hence,  impractical. 

Other  researchers  have  developed  methods  to  incorporate  explicit  geometric  information 
about  objects.  The  method  of  projection  onto  convex  sets  (POCS)  [15,16]  incorporates  prior 
knowledge  by  sequentially  projecting  candidate  estimates  onto  a  collection  of  convex  sets, 
where  each  set  represents  some  prior  knowledge.  Noise  in  the  data  tends  to  cause  the 
method  to  diverge,  however,  and  even  though  it  is  possible  to  account  for  the  noise  using 
a  smoothing  operator  [17],  finite  pixel  error  caused  by  iteration  between  Radon  space  and 
object  space  may  still  cause  convergence  to  the  wrong  solution  [14].  Other  investigators 
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avoid  the  latter  problem  by  iterating  entirely  in  projection  space  [18,19].  Another  approeich 
to  geometric  modeling  is  to  parameterize  the  object  directly  in  the  class  of  interest,  reducing 
the  number  of  p2irameters  that  must  be  estimated.  The  work  by  Rossi  and  Willsky  [20,21], 
Bresler  and  Macovski  [22,23,24],  Hanson  [25],  and  more  recently  Soumekh  [26],  Horn  [27], 
and  Fishbum  et.  al.  [28]  are  examples  of  this  kind  of  modeling. 

Our  approach  is  to  treat  the  noise,  the  physical  imaging  geometry,  and  prior  probabilis¬ 
tic  information  as  fundamenteil  and  explicitly  modeled  pieces  of  an  overzdl  inverse  problem 
formulation.  Our  solution  satisfies  the  m2iximum  a  posteriori  (MAP)  criterion  and  incor¬ 
porates  the  following  prior  geometric  information: 

•  The  values  of  line  integrals  t2iken  over  lines  close  in  either  latereil  displacement  (with 
the  same  angle)  or  angle  (with  the  S£une  lateral  displacement)  tend  to  be  similar  in 
veilue. 

•  The  Radon  transform  of  any  cross  section  obeys  certain  fundamental  mathematical 
consistency  conditions.  In  peirticular,  these  conditions  prescribe  constraints  on  the 
mass  and  center  of  mass  of  each  projection. 

•  The  convex  support  of  the  cross-section  density  function  uniquely  specifies  a  related 
region  of  support  of  the  Radon  transform. 

Since  both  the  primary  processing  and  the  introduction  of  prior  knowledge  take  pl£u:e  in 
Radon  space,  our  approach  is  a  projection-space  method.  This  takes  advantage  of  the  fact 
that  the  noise  is  well-modeled  as  white  in  this  domain,  so  that  the  criterion  of  statistical 
optimality  is  easily  specified,  but  has  the  disadvantage  that  prior  information  about  the 
object  is  not  conveniently  incorporated  in  Radon  space.  For  example,  a  local  prior  proba¬ 
bilistic  model  of  the  object  is  decidedly  non-local  when  transformed  into  Radon  space.  We 
circumvent  this  problem  by  modeling  directly  in  projection  space,  using  a  Markov  random 
field  (MRF)  model  of  sinograms  that  incorporates  the  three  properties  listed  above. 

Because  of  the  particular  form  of  the  chosen  MRF,  we  are  able  to  formulate  an  ^m^dogous 
problem  on  the  sinogram  continuum  (as  opposed  to  the  usual  lattice  system),  which  leads  to 
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a  closed-form  solution  given  by  a  partial  differential  equation  (PDE)  with  constraints.  We 
solve  this  constrained  PDE  using  a  primal-dual  optimization  approach  that  solves  the  PDE 
with  assumed  Lagrange  multipliers  in  the  primal  phase  and  updates  the  Lagrange  multipliers 
in  the  dual  phase.  The  primal  phase  is  fast  (although  iterative)  and  parallelizable,  due  to 
the  local  interactions  of  the  MRF;  the  dual  phase  is  fast  and  partially  parallelizable,  since 
it  uses  a  simple  update  formula  on  each  of  the  columns  of  the  sinogram  separately. 

The  paper  is  organized  as  follows.  In  Section  2  we  present  additional  background  related 
to  the  support  and  consistency  of  the  2-D  Radon  transform.  Section  3  develops  a  Markov 
random  field  model  of  sinograms,  and  formally  defines  the  maximum  a  posteriori  solution. 
In  Section  4,  we  present  a  fast  iterative  solution  that  solves  this  large-scale  optimization 
problem  and  Section  5  presents  some  experimental  results.  Fineilly,  we  discuss  these  results 
and  some  outstanding  problems  in  Section  6. 

2  CONSISTENCY  AND  CONVEX  SUPPORT 

2.1  Consistency  of  the  Radon  Transform 

An  important  fact  that  we  exploit  in  this  paper  is  that  not  aU  functions  p  :  IR^  x  5^  — ♦  IR^ 
are  valid  2-D  Radon  transforms.  A  valid  Radon  tr2Lnsform,  that  is,  a  function  that  is  the 
Radon  transform  of  some  function  /  :  IR*  -+  IR^,  is  constrained  to  lie  in  a  particular 
functional  subspace  of  the  space  of  all  real  functions.  This  subspace  is  characterized  by  the 
property  that  g  is  even  in  t  and  w  and  by  the  property  that  certain  generalized  Fourier 
coefficients  of  g  must  be  zero.  The  precise  mathematical  conditions  for  the  consistency  are 
given  by  the  foUowing  theorem  due  to  Ludwig  [31]. 

Theorem  1  (2-D  Consistency  Theorem)  In  order  for  g{t,  0)  to  be  the  2-D  Radon  trans¬ 
form  of  a  function  f  E  S  (IR^),  where  S  is  the  space  of  rapidly  decreasing  C°°  functions,  it 
is  necessary  and  sufficient  that 


(a)  geS{m}x  S^), 

(b)  git,9  +  it)  =  g{-t,9),  and 
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(c)  the  integral 


f  g{t,9)t^dt 

J  — OO 

be  a  homogeneous  polynomial  of  degree  k  in  cos9  and  sin 5  for  all  k>0. 
Proof  See  fSlJ,  [SO],  or  [SSj. 


(2) 

□ 


The  two  lowest  order  moments  of  g{t,  9)  give  the  mass  and  center  of  mass  constraints. 
The  mass  constraint  tells  us  that  the  integral  of  any  projection,  which  may  be  thought  of 
as  the  mass  of  the  projection,  must  have  the  same  value  for  any  9,  and  that  value  is  equal 
to  the  integral  of  f(x).  If,  for  example,  a  noisy  meeisurement  of  a  true  Radon  transform  has 
any  two  projections  that  do  not  integrate  to  the  same  value  then  the  measurement  is  not  a 
valid  Radon  transform,  and  it  follows  that  an  inverse  transform  is  theoretically  undefined. 
The  center  of  mass  constr^dnt  tells  us  that  the  (l-D)  center  of  mass  of  a  given  projection  is 
equal  to  the  projection  of  the  (2-D)  center  of  mass  of  the  object  onto  the  w-axis.  From  this 
one  can  see  that  the  collection  of  centers  of  mass  of  the  projections  for  different  9  must  be 
a  cosinusoidal  function  with  period  2^.  If  that  is  not  true  for  a  given  measurement  then, 
again,  the  measurement  is  not  a  valid  Radon  transform.  These  two  facts  are  easily  shown 
using  the  Consistency  Theorem,  and  may  be  stated  as 


m 


f  g{t,  9)dt=  f  f{x)  dx 
J-oo  Jx€R^ 


(3) 


c{9)  =  —  /  g{t,9)tdt  =  aco89 +  bsin9 ,  (4) 

TT%  J  — -00 

for  some  real  constants  o  and  b.  We  refer  to  m  eis  the  mass  of  f{x)  and  (3)  as  the  mass 
constraintfoT  the  2-D  Radon  transform;  the  quantity  c(^)  is  the  center  of  mass  of  projection 
g{t,9),  and  equation  (4)  is  the  center  of  mass  constraint.  It  is  also  true  that  the  center  of 
mass  of  the  projection  g{t,9)  is  the  projection  of  the  2-D  center  of  mass  of  /(x)  (see  [32]) 
eind,  indeed,  if  (R,  denotes  the  poleu:  coordinates  of  the  center  of  mass  of  the  object,  then 
it  can  be  shown  that  [20] 

c{9)  =  R  cos(0  —  f>)  • 
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Given  the  above  development,  we  see  that  if  the  mass  and  center  of  mass  were  known  o 
priori  then  (3)  and  (4)  should  be  imposed  as  constraints  on  the  estimated  sinogram.  The 
center  of  mass  constraint  c(0)  =  0,  imposed  in  Section  3,  implies  that  the  object  is  centered 
at  the  origin.  Given  a  known  center  of  mass  it  is  possible  to  adjust  the  observed  sinogram 
(by  shifting  each  of  the  projections  in  t)  to  make  it  appear  as  if  the  object  were  centered  at 
the  origin.  This  adjustment  may  always  be  accomplished  provided  one  has  a  field  of  view 
large  enough  to  encompass  both  the  original  object  and  the  object  shifted  to  the  origin.  We 
assume  this  to  be  the  c^lse. 

2.2  Object  Support  and  Radon  Transform  Support 

The  convex  support  of  an  object  is  smallest  convex  set  that  supports  the  function  /(x). 
In  this  section  we  develop  a  relationship  between  the  convex  support  of  em  object  and  a 
particular  region  of  support  of  its  2-D  Radon  transform.  This  relationship  is  a  special  case 
of  the  support  theorem  stated  and  proved  by  Lax  and  Phillips  in  [29]  and  also  discussed  by 
Helgason  [30]. 

Suppose  /  is  zero  outside  Dr,  the  disk  of  radius  T  centered  at  the  origin.  Then  it  is  easy 
to  see  from  the  definition  of  the  2-D  Radon  transform  in  (1)  that  g(t,  9)  must  be  zero  when 
t  ^  [-r,T].  Using  the  periodicity  of  the  2-D  Radon  transform  established  in  Theorem  1, 
one  can  now  conclude  that  g{t,  0)  is  completely  determined  by  its  values  on  the  rectangle 

Vt  =  {it,0)  \  -T<t<T,0<9<ir}.  (5) 

But  this  idea  can  be  refined  even  further.  Let  T  be  the  set  of  points  in  Yp  for  which  /(x)  #  0. 
Now  consider  the  Radon  transform  g{t,  9)  of  /,  and  the  unit  vector  oj  =  [cos^  sin^]"^.  With 
reference  to  Fig.  2  and  to  (1),  we  see  that  for  any  given  w,  the  value  of  the  Radon  transform 
must  be  zero  for  t  >  2ind  t  <  .  Here,  is  the  lateral  position  of  the  line  perpendicular 

to  u  which  is  positioned  as  far  as  possible  in  the  -f-w  direction  so  it  just  grazes  the  set  T ; 

is  the  lateral  position  of  the  line  perpendicular  to  w  which  is  positioned  as  far  a  possible 
in  the  -w  direction  so  it  just  grazes  the  set  7.  The  quantities  i+  and  t-  are  called  support 
values  amd  the  corresponding  lines  are  called  support  lines  of  the  set  7.  Knowledge  of  both 
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t+  and  t-  for  all-in  [0,  ?r)  determines  the  convex  hull  of  7,  denoted  hul(7),  which  is,  by 
definition,  the  smallest  convex  set  cont2Lining  7.  The  set  hul(7)  is  also  the  convex  support 
of  the  function  /(*). 

From  the  above  discussion,  we  conclude  that  the  2-D  Radon  transform  is  completely 
determined  by  its  values  on  the  set 

g  =  {{t,9)eyT\t^{e)<t<t+{9)}  (6) 

where,  for  clarity,  we  have  explicitly  indicated  the  functional  dependence  of  t+  and  t-  on 
9.  An  example  of  such  a  set  is  shown  in  Fig.  3.  For  a  given  object  support  set  7,  we  think 
of  ^  as  the  matching  region  of  support  in  Radon  space.  However,  although  7  uniquely 
determines  it  is  clear  that  Q  determines  only  hul(7),  not  7  itself.  Furthermore,  g  is 
not  necessarily  the  actual  support  of  g(t,9)  since  it  is  possible  for  g{t,9)  to  be  zero  when 
(t,9)  G  g  if  T  is  not  connected.  We  eure  primarily  concerned  with  the  convex  support  of  /, 
since  this  is  what  may  be  determined  directly  from  knowledge  of  g. 

In  Section  3  we  assume  that  an  estimate  of  g  is  available,  and  we  define  a  prior  prob¬ 
ability  on  sinograms  that  gives  a  low  probability  to  sinograms  that  have  non-zero  values 
outside  of  g. 

3  SINOGRAM  MRF  AND  MAP  ESTIMATION 

We  chose  to  represent  prior  probabilistic  knowledge  about  sinograms  using  a  Markov  random 
field  (MRF)  on  a  discrete  sinogram  lattice.  There  are  several  reasons  for  this  choice.  First, 
the  MRF  is  a  convenient  way  to  describe  processes  with  local  interaction  in  such  a  way  that 
the  joint  probability  over  all  sites  is  easily  determined.  Second,  the  constraints  that  arise 
from  the  2-D  Radon  transform  consistency  conditions  are  easily  incorporated  in  the  MRF 
by  limiting  the  space  of  allowable  configurations.  Third,  knowledge  of  the  convex  support 
of  the  object,  which  is  treated  £is  a  penalty  rather  than  a  constreunt,  may  be  incorporated 
into  the  MRF  by  adding  an  additional  self-potential  term  (see  below).  FineJly,  this  choice, 
along  with  the  particular  details  described  below,  leads  directly  to  a  statistically  optimum 
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maximum  a  posteriori  solution.  We  will  see  in  the  following  section  that  the  form  of  this 
estimate  has  an  analogous  variational  formulation  that,  once  solved,  leads  to  a  fast  iterative 
solution. 


3.1  A  Sinogram  MRF 

This  section  develops  a  Markov  R^lndom  field  (MRF)  on  the  sinogram  lattice  that  includes 
the  mass  and  center  of  mass  constraints,  and  that  includes  the  periodicity  and  smoothness 
of  the  2-D  Radon  transform  and  the  convex  support  of  the  object.  The  ingredients  needed 
to  define  a  MRF  are  [33]:  1)  the  lattice,  2)  the  potential  functions,  3)  the  graph  structure, 
eind  4}  the  feasible  configurations. 


The  Sinogram  Lattice  As  discussed  in  Sections  1  and  2,  a  sinogram  is  an  image  of  the 
Radon  transform  (or  meEisured  Radon  transform)  of  an  object  over  the  truncated  domain 
1/r  (see  equation  (5)),  with  brighter  intensities  corresponding  to  larger  values  of  g{t,6). 
In  order  to  define  a  MRF,  however,  we  require  a  finite  lattice  system,  rather  than  the 
continuum  of  points  in  l/j.  Therefore,  we  define  the  sinogram  lattice  to  be  a  rectangularly 
sampled  version  of  j/r  given  by 


i/5  =  {(*,<>)  I  < 


n<i  -  1 
2 


—  1 
2  ’ 


0, 


(7) 


where  nj  is  the  number  of  sample  points  in  t,  and  n„  is  the  number  of  sample  points  in  9. 

For  convenience  we  adopt  the  following  notation  for  the  remainder  of  this  section.  A 
site  in  the  sinogram  lattice  is  denoted  s  =  (i,j)  and  the  set  of  all  sites  by  5.  A  site  value 
is  denoted  in  several  ways:  g,  —  gij  —  g{ti,9j).  The  collection  of  all  site  values  is  called 
the  discrete  sinogram  or  just  the  sinogram  when  the  meaning  is  cleu  from  the  context. 
Note,  that  a  site  in  the  sinogreun  lattice  corresponds  to  a  line  in  the  pl2ine  passing  through 
the  disk  Dt-  In  particular,  the  site  (i,y)  corresponds  to  the  line  L{ti,9j)  —  {(a:, y)  € 
IR*  I  xcoa9j  +  yam9j  =  t,}. 


The  Potential  Functions  The  physics  of  this  problem  does  not  specify  for  us  a  prior 
probability  on  sinograms.  We  rely  on  experimentation  and  intuition  to  surmise  what  a 
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reasonable  form  for  a  prior  might  be,  given  what  we  know  about  the  types  of  objects  under 
consideration  and  the  transformation  of  those  objects  via  the  Radon  transform.  After  much 
thought  and  experimentation,  one  key  idea  has  driven  us  to  implement  what  turns  out  to  be 
the  simplest  kind  of  MRF.  This  idea  is  simply  that  sinograms  tend  to  be  locally  smooth.  A 
prior  that  produces  such  sample  functions  is  an  MRF  with  potential  functions  that  prescribe 
an  affinity  between  nearest  neighbors  —  this  is  the  so-called  nearest-neighbor  “blob”  model 
[34].  Let  s,r  €  S  denote  sites  that  are  either  vertical  or  horizontal  nearest  neighbors.  To 
prescribe  affinities  between  the  sinogram  values  defined  on  these  sites  we  define  the  vertical 
pair-potential  function  as 

—  ^v{9$  ~  <7r)*  (8) 

and  the  horizonted  pair-potential  function  as 

—  ^h{9*  ~  9r)^  (9) 


where  (s,  r)  and  {s,  r)  represent  pairs  of  adjacent  sites  in  the  vertical  and  horizontal  direc¬ 
tions,  respectively.  The  positive  constants  b„  and  6/,  allow  one  to  make  the  vertical  and 
horizontal  affinities  of  different  strength,  thus  making  this  a  non-isotropic  random  field  [35]. 
In  this  paper,  we  choose  the  constants  b„  and  bf^  a  priori]  however,  it  is  possible  to  use  the 
actual  data  to  estimate  these  coefficients  as  part  of  a  hierarchical  estimation  algorithm  [36]. 

The  self-potentials  are  defined  using  knowledge  of  the  object’s  support.  As  developed  in 
Section  2,  the  object’s  support  7  implies  a  matching  region  of  support  Q  within  the  Radon 
transform  domain  t/y.  If  either  set  were  known  exactly  then  we  would  insist  that  sinogram 
values  in  the  region 

9  =  yT-9  (10) 


be  exactly  zero.  However,  we  shall  assume  that  only  an  estimate  of  the  sinogram  support 
is  available,  and  that  we  have  a  measure  of  that  estimate’s  reliability.  Therefore,  sinograms 
with  non-zero  values  in  9  should  have  low  probability,  but  not  as  low  if  the  support  estimate 
is  unreliable.  To  provide  this  effect,  we  define  the  support  self-potential  as 


otherwise 


(11) 
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where  gg  stands  for  the  value  of  the  sinogram  at  site  s=(i,j),  and  k  is  a  positive  constant 
which  is  used  to  reflect  the  support  measurement’s  reliability. 

The  Graph  Structure  The  form  of  the  potential  functions  described  above  dictate  the 
required  neighborhood  structure,  and,  in  feict,  only  nearest  neighbors  aire  necessary.  In  this 
case,  the  most  general  form  of  the  MRF  energy  function  is  [33] 

U{9)  =  H  y{»,r){9»,9r)  +  X)  (12) 

*  (».'•)  <*.»•> 

where  the  first  summation  is  taken  over  all  sites  in  the  lattice,  the  second  over  all  vertical 
nearest-neighbors,  and  the  third  over  aJl  horizontal  nearest-neighbors. 

Since  all  objects  are  known  to  be  zero  outside  the  disk  of  radius  T,  the  boundary  value 
above  and  below  the  sinogram  must  be  zero  as  shown  in  Fig.  4a.  The  boundaries  to  the  left 
and  right  of  the  sinogram  must  be  treated  difierently,  however.  Here,  we  use  the  symmetry 
property  of  the  consistency  theorem  stated  in  Section  2 

9{t>9)  =  g{-t,0  +  v) 

to  conclude  that  the  neighbors  wrap  around  in  a  toroid  that  is  twisted  or  flipped  about  the 
t-axis  as  shown  in  Fig.  4b.  In  other  words,  the  sinogram  is  actually  defined  on  a  Mobius 
strip.  Given  these  boundaries  conditions,  peiir  potentials  involving  one  site  outside  of  the 
lattice  may  now  be  computed.  Sites  above  or  below  the  sinogram  are  given  the  fixed  value 
zero,  whereas  sites  to  the  left  or  right  are  given  the  value  of  the  sinogram  on  the  opposite 
side  of  the  sinogram  and  in  the  opposite  t  direction. 

The  Feasible  Configturations  The  mass  and  center  of  mass  constraints  used  for  the 
MRF  are  discrete  approximations  to  the  integral  expressions  of  (3)  and  (4).  Thus,  letting 
m  denote  the  object’s  mass  we  have 

^  =  ">  Vj,  1  <  y  <  n„  (13) 

and 

=  0  Vj,  l<j<n„  (14) 

Si 
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where 


2T..  nrf  +  1 
U  =  — — ir~ 

n<i  I 

is  the  lateral  position  of  the  tth  line. 

The  presence  of  constraints,  even  linear  equality  constraints  such  eis  these,  makes  the 
computation  of  the  MAP  estimate  more  di£Scult  since  the  algorithm  must  be  a  constrained 
optimization  method  [37].  In  fact,  the  solution  must  be  an  element  of  a  set  of  feasible 
discrete  sinograms  Oj,  which  contains  all  real  matrices  of  dimension  by  that  satisfy 
(13)  and  (14). 


3.2  The  MAP  Formulation 

The  Gibbs  Prior  Having  now  defined  all  the  elements  of  the  MRF,  the  joint  probability 
density  for  the  discrete  sinogram  prior  is  simply  given  by  the  Gibbs  density 


p{9)  =  5  €  n,  (16) 

where  g  denotes  the  vector  of  sinogram  site  values  and  Z  is  given  by 


so  that  p{g)  integrates  to  one.  The  function  U{g)  is  the  energy  function  defined  in  (12). 


The  Observations  We  assume  that  noisy  observations  of  the  true  site  values  are  available 
over  a  (possibly)  limited-angle  or  sparse-angle  subset  t/o  of  yx  and  that  the  observations 
are  given  by 

y»y  =  9ij  +  (t«  j  ^y)  ^  Vo  >  (l^) 

where  the  rnj  are  independent  zero-mean  Gaussian  rEmdom  variables  with  variance  <7^. 
Letting  g  denote  the  vector  of  true  sinogram  site  values  and  letting  y  and  n  denote  the  vector 
of  observations  and  noise  samples,  respectively,  we  may  write  the  observation  equation  in 
vector  form  as 

y  =  Sg->rn  (18) 
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where  S  is  a  matrix  that  selects  the  observations  as  foUows.  In  the  measurement  configu¬ 
ration  we  consider,  a  column  of  the  matrix  given  by  [y,-,-]  is  either  observed  completely  (in 
additive  noise)  or  not  observed  at  all.  Suppose,  for  the  purposes  of  this  discussion  only,  we 
form  g  by  stacking  the  colunms  of  [j,-,],  stacking  2J1  the  observed  columns  first,  from  the  top 
proceeding  downwards,  and  the  remaining  columns  following  in  any  order.  Then,  denoting 
the  number  of  observed  columns  by  n^,  we  see  that  5  is  given  by 

s  =  [i\o] 

where  I  is  the  no^d  x  riorid  identity  matrix. 

The  Sinogram  MAP  Estimate  Now,  with  the  observation  equation  given  by  (18)  and 
the  prior  probability  given  by  (16),  we  may  now  derive  the  form  of  the  MAP  estimate 
9map-  Denoting  the  noise  covariance  matrix  by  Kn  we  may  use  (18)  to  write  the  conditional 
measurement  density  (zero  mean,  jointly  Gaussian)  as 

P(y\9)  =  |2z- exp  ^(y  -  S9)’^K-^{y  -  Sg)^  .  (19) 

Using  the  definition  of  conditional  probability  twice,  and  the  prior  probability  given  by  (16), 
the  posterior  distribution  is  found  to  be 

p(y|y)  =  p{yiy)p(5)/p(y)  (20) 

=  gxp(-(^(y  -  -  ^9)  +  u{g))) 

The  MAP  estimate,  which  maximizes  (20)  with  the  true  observations  Y  substituted  into 
the  expression,  is  given  by 

ymap  =  argmin^(F-5j)'^(F-55)-l-17(s)  (21) 

gertf  2ff' 

where  we  have  used  the  fact  that  Kn  =  cr^  I. 

The  posterior  distribution  of  (20)  is  a  Gibbs  density,  and  since  the  observation  equation 
is  not  convolutional,  its  graph  structure  is  identical  to  that  of  the  prior.  The  only  significant 
difiference  between  the  two  MRF’s  is  the  form  of  the  energy  function,  which,  in  the  posterior 
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density,  contains  a  self  potential  term  that  couples  the  observations  y  to  the  sinogram  g. 
The  identification  of  the  posterior  density  as  a  Gibbs  density  serves  as  the  basis  for  the 
stochastic  relaxation  and  simulated  annea/tny  methods  of  Geman  and  Gemein  [33]  and  others, 
eilgorithms  that  have  been  the  focus  of  much  research  in  recent  years.  These  methods, 
besides  being  generally  very  slow,  are  inappropriate  for  this  application  for  two  reasons. 
First,  the  stochastic  methods  do  not  conveniently  incorporate  construnts  [32].  Second, 
and  most  importantly,  the  minimization  problem  of  (21)  requires  minimizing  a  quadratic 
function  with  linear  constraints,  which  when  taken  advantage  of  as  we  do  in  the  next  section, 
leads  to  a  much  faster  algorithm. 

4  Primal-Dual  MAP  Algorithm 

This  section  develops  the  theory  and  implementation  of  a  fast  iterative  algorithm  for  com¬ 
puting  gmap-  The  key  step  in  the  development  this  method  is  to  write  the  vector  minimiza¬ 
tion  problem  of  (21)  as  a  construned  minimization  problem  involving  an  unknown  function 
g{t,  9)  over  the  continuous  domain  I/t.  The  solution  to  this  variational  formulation  is  a 
partial  differential  equation  (PDE)  with  three  unknown  functions:  the  sinogram  and  two 
Lagrange  multiplier  functions.  We  use  a  generic  primal-dual  method  to  find  the  solution  to 
this  PDE,  incorporating  a  fast  iterative  local  relaxation  algorithm  in  the  primal  stage  and 
simple  Lagrange  multiplier  update  formulas  in  the  dual  stage. 

4.1  Variational  MAP  Formulation  and  Solution 

The  Miulmizatiozi  Problem  Consider  the  problem,  which  we  refer  to  as  (V),  of  mini¬ 
mizing 

‘ = I L  If}  IL  K^)  ^  (If )  ]  * 

subject  to  the  equality  constreiints 

Jl  =  m  =  /  g{t,  9)  dt  (23) 

J-T 

1 

J2  =0=  —  /  tg{t,9)dt 

m  J-T 
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and  boundary  conditions 


5(T,^>)  =  9{-T,0)  =  O 
9it,0)  =  9i-t,^) 


(24) 


where  k,  and  7  are  positive  constants.  This  problem  is  a  continuous  formulation  of  the 
sinogram  MAP  problem  of  (21).  The  first  term  in  I  is  analogous  to  the  first  term  in  (21) 
—  both  represent  a  penalty  that  seeks  to  keep  the  estimate  close  to  the  observations.  The 
second  two  terms  are  analogous  to  the  terms  of  17(^),  given  in  (12).  The  first  term  comprises 
the  support  information,  which  matches  the  summation  over  single  sites  in  U{g).  The  second 
term  integrates  the  sum  of  the  squares  of  the  two  partial  derivatives  of  9,  which  corresponds 
to  the  two  summations  of  peur-potentials  in  U{9).  The  two  integral  constraints  in  (23)  are 
exactly  the  mass  and  center  of  mass  constraints.  Finally,  the  boundary  conditions,  which 
include  the  twisted  boundary,  we  stated  in  (24). 

To  simplify  notation  we  define  the  following  indicator  function  notations: 


XG{t,e)  = 


1  {t,9)eg 

0  otherwise 

which  indicates  the  complement  of  the  region  of  support  in  the  sinogram  domain,  and 


(25) 


: 


(*,«)€  I/O 
otherwise 


(26) 


which  indicates  the  region  in  the  sinogram  over  which  observations  are  available.  Using  this 
notation  we  may  write  I  as 


The  problem  is  now  in  the  form  of  a  classical  variational  problem  which  may  be  solved  using 
standard  calculus  of  variations  techniques  (see  [38],  for  example). 


Partial  DifiTerential  Equation  A  necessary  condition  for  9{t,0)  to  be  a  solution  to  (V) 
is  that  it  satisfy  the  following  second  order  partial  differential  equation  (PDE)  [32] 

(2kXg  +  jjXr)  9  -  2^0  -  27!^  =  ^XyV  -  Ai(5)  -  X2i0)t  (28) 
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and  the  additional  boundary  condition 


d9{t,0)  ^  dg{-t,ic) 
dt  dt 


In  addition,  g{t,  9)  must  satisfy  the  original  constrmnts  and  boundary  conditions.  Note  that 
(28)  contedns  three  unknown  functions:  g{t,  0),  and  two  Lagrange  multiplier  functions  Ai(0) 
and  A2(^)  (one  for  each  constraint).  To  simplify  the  expressions  in  the  remainder  of  this 
section  we  use  the  notation  gt  and  ga  to  stand  for  the  first  and  second  partial  derivatives 
of  g{t,9)  with  respect  to  t,  respectively,  and  9$  and  9$$  to  stand  for  the  first  and  second 
partial  derivatives  of  g{t,0)  with  respect  to  9,  respectively. 

To  solve  (28)  for  g{t,9),  Ai  and  A2  must  first  be  determined.  An  analytic  expression 
for  Xi  may  be  found  by  integrating  both  sides  of  (28)  zmd  simplifying;  and  an  analogous 
expression  for  A2  may  be  found  by  multiplying  both  sides  of  (28)  by  t,  then  integrating  and 
simplifying.  The  results  are  [32] 


Ai(0) 

A2(<>) 


2T 

-3 

2T3 


f  2kXg9  dt  -  2pgt\^T  + 
J-T 

rT 

I  2tKXG9  dt  -  2^tgt\_j< 
J-T 


m  1  ^ 

1 

— 2  /  ^^rydt  , 
flT-*  J-T 


(30) 

(31) 


These  equations  may  be  substituted  into  (28)  to  give  an  integro-differential  equation  in  a 
single  unknown  function  g{t,9). 


Primal-Dual  Optimization  Method  Unfortunately,  although  the  resulting  integro- 
difierential  equation  could  be  discretized  and  solved  numericedly,  the  problem  is  very  leirge 
emd  computationally  intractable.  The  approach  we  take  instead  is  to  use  the  generic  primal- 
dual  method  described  by  Bertsekas  in  [39]  to  solve  the  PDE  directly.  In  outline,  the 
method  requires  us  to  solve  (28)  numerically  given  estimated  Lagrange  multipliers,  update 
the  Lagrange  multipliers  if  the  solution  doesn’t  meet  the  required  constraints,  and  repeat 

A  A 

until  a  jointly  optimum  trio  of  g,  Ai,  and  A2  is  found. 

To  find  an  initial  estimate  of  the  Lagrange  multipliers  we  make  several  approximations 
which  often  hold  true  at  or  near  the  solution.  First,  near  the  solution  we  expect  that 
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g{t,  w  0  for  (t,  9)  e  especially  when  k  is  large.  Hence,  we  may  make  the  approximations 

f  2kXg9  dt  mO  and  f  2tKXG9dt  w  0. 

J-T  J-T 


Second,  the  terms  in  (30)  amd  (31)  involving  the  partial  derivative  9t  evaluated  at  ±T  may 
often  be  close  to  be  zero.  These  approximations,  applied  to  (30)  and  (31),  yield  the  following 
initial  Lagrange  multiplier  estimates 

^  (32) 

^  ’  (23) 

each  of  which  may  be  evaluated  given  only  the  data.  Substituting  these  functions  for  the 
true  Lagrange  multipliers  in  (28)  yields  the  PDE 


^2kXg:  +  9  -  ^P9u  -  2'f9e0  = 

1  1  /-^  j  m  Zt 

^  L  ~  J-T  ’ 


(34) 


which,  unlike  the  original,  has  a  single  unknown  9{t,0),  and  may  be  solved  numerically 
using  any  of  several  techniques  for  solving  elliptic  PDEs  as  discussed  below. 

If  g  does  not  meet  the  constraints  after  solving  (34),  the  first  primal  stage,  then  we 
conclude  that  the  approximations  used  to  derive  the  approximate  PDE  of  (34)  were  not 
accurate.  This  situation  will  in  fact  occur  if  ^  is  leirge  or  if  the  observations  are  missing 
entire  projections,  such  as  in  the  limited-angle  and  sparse-angle  problems.  Therefore,  the 
dual  stage  updates  the  Lagrange  multiplier  functions  to  move  them  closer  to  their  optimum 
values  using  the  following  formulas  [39] 

=  A{(«l)+aU-£^a(i,»)<«j  (35) 

A*+-(«)  =  A‘(«)+a(’o-i|_%,(t,9)*l  (36) 

where  a  is  a  positive  constant,  and  k  is  an  iteration  counter.  Note  that  the  update  cal¬ 
culations  can  be  done  independently,  and  therefore  in  parallel,  for  each  projection  in  the 
sinogram.  After  each  update,  the  new  Lagrange  multiplier  functions  are  substituted  into 


(28),  which  is  solved  numerically  for  a  new  g.  When  g  meets  the  constraints  to  within  a 
specified  tolerance,  then  the  three  functions  are  jointly  optimal  and  g  is  the  desired  sinogram 
estimate. 

The  constant  a,  which  appears  in  (35)  and  (36),  is  chosen  l^u:ge  enough  so  that  conver¬ 
gence  to  the  correct  Lagrange  multipliers  (and,  hence,  the  correct  solution  to  (V))  occurs 
quickly,  yet  not  so  large  that  the  sequence  will  not  converge.  Bertsekas  [39]  describes  the 
the  selection  of  a  and  relates  this  generic  primal-dual  method  to  the  method  of  multipliers, 
about  which  a  great  desil  of  theory  is  known.  In  our  experiments,  a  was  chosen  empirically 
to  yield  a  good  rate  of  convergence  for  our  problem. 


4.2  Numerical  Methods 

The  sinogram  is  approximated  on  a  rectilinear  grid  with  vertical  and  horizonted  sample 
spacings  given  by  At  =  2Tfnd  and  Ag  =  Jr/fiu,  respectively.  The  PDE  of  (28)  may  then  be 
approximated  at  an  interior  point  by  the  finite  difiference  equation  [40] 


where 


U,j9i-\j  ~  U,i9ij+i  —  bijgij~i  —  Sij 


(37) 


Uj  =  2^  (38) 

ri.i  =  2$ 


Ki  =  27 

ti,i  =  27 

dij  =  4(9  -f  47  -1-  (2kXg  +  I 

and  /9  =  At  and  7  =  7/Aj.  Sinogreun  values  in  (37)  that  correspond  to  points  outside 
the  lattice  must  be  evaluated  according  to  the  boundary  conditions  developed  in  Section  3. 
The  set  of  equations  given  by  (37)  for  all  j,  j  =  1, . . . ,  nj  and  i,i  =  1, . . . ,  n„  may  be 
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organized  and  written  as  a  vector  equation 

Ag  =  s  ,  (39) 

and,  although  this  is  a  very  high  dimensional  problem,  the  local  interactions  that  result  from 
the  nearest-neighbor  construction  allow  for  efficient  iterative  solutions.  Several  traditional 
methods  (cf.  [41])  including  Jacobi,  simultaneous  over-relaxation  (SOR),  and  Chebyshev 
semi-iterative  relaxation  methods  may  be  employed  to  solve  (39).  However,  we  have  chosen 
to  implement  a  relatively  new  method  due  to  Kuo,  Levy,  and  Musicus  [40]  which  has  been 
shown  to  have  very  favorable  convergence  properties,  and  is  relatively  easy  to  implement. 
This  method,  in  addition,  has  been  shown  to  be  ideally  suited  for  parallel  implementation. 

Our  implementation  of  Kuo’s  local  relaxation  algorithm  (KLR)  is  a  special  case  of  the 
more  general  implementation  described  in  [40] .  We  assume  the  PDE  to  be  of  the  form 

d^u  d^u  ,  f  ^  ff  ^ 

~ 

where  (®i,X2)  G  [0, 1]  x  [0, 1],  and  to  satisfy  the  conditions  given  in  [40].  Then  the  PDE  is 
approximated  by  the  5-point  stencil 


with 


l  =  P,  r  =  p,  b  =  q,  t  =  q, 
di,}  =2p  +  2q  +  Sij  =  fijh'^ 


where  h  is  the  grid  spacing  and  is  defined  as  i{ih,jh).  Each  grid  point  is  assigned  a 
color,  either  red  or  black,  2u:cording  to  an  alternating  pattern  as  on  a  checkerboard.  Then 
the  local  relaxation  procedure  can  be  written  as: 


red  points  (i  +j  is  even): 

=  (1  -  +  «<.y)  . 

black  points  (*  +j  is  odd): 

=  (1  -  +  Uijdrj  > 
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where  is  called  the  local  relaxation  parameter  and  is  given  by 


2 


where 

2  (  ;r  .  JT  \ 

One  point  related  to  convergence  of  KLR  is  worthy  of  comment.  In  the  initialization 
phase,  KLR  calculates  an  array  of  local  relaxation  parameters,  coij,  one  per  site,  which 
are  theoretically  optimum  for  a  particular  boundary  condition  which  our  problem  does 
not  satisfy  (because  of  the  twisted  boundary  property).  Therefore,  KLR  still  converges  to 
the  correct  solution,  but  it  may  do  so  more  slowly  than  the  predicted  convergence  rates. 
However,  we  have  found  in  our  experiments  that  no  slow-down  is  evident;  the  rate,  in 
practice,  is  still  of  order  y/N,  where  N  is  the  total  number  of  points  in  the  grid. 


5  Experimental  Results 

5.1  Overview 

In  this  section,  we  present  results  from  several  simulation  studies,  each  designed  to  demon¬ 
strate  a  different  eispect  of  the  sinogram  MAP  algorithms  described  in  Section  4.  The  object 
that  is  used  for  all  of  the  simulations  in  this  section  is  an  ellipse  with  the  letters  M  I  T  in 
its  interior,  as  shown  in  Fig.  5.  The  ellipse  is  centered  at  the  origin  and  rotated  45  degrees 
in  the  clockwise  direction,  and  has  two  values:  0  outside  of  the  ellipse  and  1  within  the 
body  of  the  ellipse,  except  within  the  letters,  where  the  vzilue  is  0.  The  noise-free  sinogram 
shown  in  Fig.  6a  is  calculated  using  approximate  strip  integrals  (see  [42])  from  analytic 
expressions  for  the  ellipse  and  characters  in  the  interior.  The  sinogram  has  81  rows  and 
60  columns,  approximating  g{t,$)  over  the  angular  range  [x/2,3jr/2)]  and  the  lateral  range 
[— r,T].  The  81x81  pixel  image  in  Fig.  6b  is  a  reconstruction  from  the  noise-free  data  of 
Fig.  6a  using  convolution  backprojection  (CBP)  with  a  ramp  filter  (see  [42]). 

The  MIT  ellipse  was  chosen  as  a  test  object  for  this  experiment  because  the  loss  of 
angular  information  has  strikingly  different  effects  depending  on  where  the  missing  angles 
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occur.  For  example,  if  the  missing  projections  have  lines  that  integrate  along  the  long 
axis  of  the  ellipse,  then  the  narrowness  of  the  ellipse  cannot  be  observed,  but  the  detail 
of  the  letters  in  the  interior  is  quite  apparent  from  the  observed  projections.  If,  however, 
the  missing  angles  occur  along  the  short  axis  of  the  ellipse,  then  the  letters  cannot  be 
observed  well,  but  the  narrowness  ta  apparent.  The  first  case  is  where  information  about 
the  boundary  of  the  object  has  a  striking  effect  on  the  reconstruction;  the  smoothing  effect 
helps  in  both  cases. 

One  noisy,  two  limited-angle,  and  two  sparse-angle  cases  were  derived  from  the  noise-free 
sinogram  and  used  as  simulated  observations.  Fig.  6c  shows  a  noisy  sinogram,  created  by 
adding  independent  samples  of  zero-mean  Gaussian  noise  with  variance  er*  to  eeich  element 
of  the  true  sinogram  of  Fig.  6a.  The  signal  to  noise  ratio  (SNR)  of  this  sinogram  is  S.OdB, 
using  the  following  definition  of  SNR: 


'  T.T.sHu.ti) 


SNR  =  10  log 


(40) 


where  g{ti,0j)  is  the  true  sinogreim.  Fig.  6d  shows  a  reconstruction  of  Fig.  6c  using  GBP. 
Figs.  7a  £ind  7b  show  speirse-angle  reconstructions  from,  respectively,  15  and  10  evenly 
spaced  projections  of  a  lO.OdB  noisy  observed  sinogram  of  the  ellipse  (not  shown).  This 
corresponds  to  projections  t2iken  12.0  and  18.0  degrees  apart,  respectively.  Figs.  7c  and  7d 
£ire  reconstructions  from  the  first  (left)  40  projections  and  the  last  (right)  40  projections, 
respectively,  of  the  same  lO.OdB  noisy  sinogreun.  The  first  limited-angle  arrangement  lacks 
projections  with  information  about  the  narrow  dimension  of  the  ellipse,  while  the  second 
arrangement  lacks  detailed  information  about  the  letters  within  the  ellipse.  Each  of  these 
four  reconstructions  was  made  using  GBP,  assigning  zero  to  the  missing  projections. 

Given  an  observed  sinogram  (perhaps  noisy  or  only  partially  observed),  the  first  pro¬ 
cessing  step  is  to  estimate  the  object  mass  using 


m 


€yrp _ 


Jni 


ye/f=i 


(41) 


where  J  is  the  set  of  observed  projections  and  J  is  the  number  of  elements  in  J.  Each 
sinogram  is  then  normalized  by  dividing  each  observation  by  m  so  that  the  normalized 
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sinogram  corresponds  to  an  object  with  unit  mass.  This  normalization  is  necessary  so  that 
the  coefficients  p,  7,  and  k  have  the  same  qualitative  effect  on  low  mass  sinograms  as  on 
high  mass  sinograms.  Using  the  normalized  sinogram,  the  noise  variance  is  estimated  from 
the  top  and  bottom  rows,  and  this  estimate  is  used  as  the  true  variance  in  subsequent 
computations. 

In  general,  the  center  of  mass  must  also  be  estimated,  perhaps  using  methods  described 
in  [32]  or  [21],  so  that  the  observed  projections  can  be  shifted  to  correspond  to  an  object 
centered  at  the  origin.  In  these  simulations,  however,  we  assume  that  the  object  is  already 
centered  at  the  origin  —  which  is  very  nearly  true  for  the  MIT  ellipse.  In  order  to  study 
the  effect  of  correct  and  incorrect  convex  hull  estimates,  the  convex  hulls  used  in  these 
studies  are  fixed  and  known,  although  not  always  correct.  Experiments  that  show  the  full 
hierarchical  procedure  that  first  estimates  the  convex  hull  are  described  in  [32] . 

5.2  Effect  of  Smoothing  Coefficients 

The  coefficient  7  has  the  effect  of  smoothing  or  blurring  the  sinogram  in  the  horizonteJ 
direction;  the  coefficient  has  a  simileur  smoothing  effect  in  the  vertical  direction.  Fig.  8 
shows  sinogram  MAP  estimates  resulting  from  the  full-view  observations  of  Fig.  6c,  using 
no  support  information.  Fig.  8a  corresponds  to  7  =  0.05  and  P  =  0.01,  Fig.  8b  to  7  =  0.5 
and  P  =  0.01,  Fig.  8c  to  7  =  0.05  and  P  =  0.1,  and  Fig.  8d  to  7  =  0.005  and  p  =  0.001. 
Images  reconstructed  from  these  sinograms  using  CBP  are  shown  in  the  respective  panels 
of  Fig.  9.  The  reconstruction  in  Fig.  9a  —  which  used  what  has  empirically  shown  to  be 
good  smoothing  coefficients  —  should  be  compared  to  the  unprocessed  CBP  reconstruction 
of  Fig.  6d. 

It  should  be  noted  from  Fig.  9b  that  excessive  smoothing  of  the  sinogram  in  the  horizon¬ 
tal  direction  results  in  circular  blurring  of  the  reconstructed  image.  Similarly,  the  haziness  of 
the  image  in  Fig.  9c  results  from  excess  smoothing  of  the  sinogram  in  the  vertical  direction, 
which  effectively  produces  a  low-pass  filtering  effect  on  each  projection.  There  is  notice¬ 
able  improvement  in  both  reconstructions  shown  in  Figs.  9a  and  9b  over  that  in  Fig.  6d; 
however,  there  are  important  differences.  For  one,  the  contrast  between  the  ellipse  body 
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and  the  background  is  better  for  the  larger  smoothing  coeflRcients  of  Fig.  9a.  However,  that 
enheincement  also  2iccompanies  a  decreased  definition  of  the  ellipse  boundary.  The  legibility 
of  the  internal  letters,  however,  appears  to  be  best  in  the  highest  contrast  image  shown  in 
Fig.  9a. 

5.3  Effect  of  Known  Support 

Fig.  10  shows  the  effect  of  varying  k  for  known  (correct)  support.  The  different  values  of 
K  are  given  by  (a)  k  =  0.0,  (b)  k  =  5.0,  (c)  k  =  10.0,  and  (d)  k  =  10,000.  In  each  case, 
the  full-view  observations  of  Fig.  6c  were  used,  and  7  =  0.05  and  jd  =  0.01.  The  object 
reconstructions  were  made  using  full- view  CBP,  and  should  be  compared  to  those  of  Figs.  6d 
and  9a,b,c,d. 

We  see  from  the  set  of  experiments  shown  in  Fig.  10  that  known  support  sharpens  the 
boundary  of  the  ellipse  considerably.  However,  in  the  image  with  the  sharpest  boundary 
(Fig.  lOd),  the  letters  in  the  ellipse  are  not  as  legible  as  the  images  in  the  other  panels  — 
the  contrast  of  the  letters  does  not  appear  to  be  as  large.  This  is  likely  to  be  due  to  the 
mass  constraint,  which,  for  k  large,  must  produce  an  estimate  that  has  all  its  mass  (for  a 
given  projection)  between  the  two  support  values.  But,  in  addition  there  is  a  smoothness 
requirement  which  is  attempting  to  reduce  abrupt  variations  within  the  projections.  This 
may  have  the  overall  effect  of  increasing  the  magnitude  of  (normally  small)  values  of  line 
integrals  that  pass  through  the  internal  letters. 

5.4  Effect  of  Incorrect  Support 

In  this  set  of  experiments  we  examine  the  effect  of  using  support  information  which  is 
incorrect.  Fig.  11  shows  results  where  the  support  corresponds  to  an  ellipse  which  has  been 
rotated  90  degrees  from  the  correct  orientation.  The  observed  sinogram  is  that  of  Fig.  6c, 
and  the  algorithm  used  the  smoothing  coefficients  7  =  0.05  and  P  =  0.01.  The  different 
reconstructions  in  Fig.  11  correspond  to  setting  k  to  (a)  k  =  0.0,  (b)  k  =  5.0,  (c)  «  =  10.0, 
and  (d)  k  =  10, 000. 

This  set  of  experiments  shows  that  as  k  grows  larger,  the  image  values  outside  the 
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assumed  region  of  support  grow  smaller.  Eventually,  this  effect  overwhelms  the  evidence 
of  the  observations  and  virtually  obliterates  the  parts  of  the  true  ellipse  that  lie  outside  of 
the  incorrect  support.  But  the  effect  of  the  mass  constraint  and  the  smoothing  coefficients 
also  affect  the  appearance  of  the  final  image.  Since  each  projection  has  mass  m,  when  the 
support  width  is  incorrectly  narrowed,  and  k  is  too  large,  then  the  sinogram  values  must 
be  very  large  within  the  region  of  support  just  to  accommodate  the  required  mass,  and  the 
values  will  typically  be  very  much  larger  than  the  observations.  As  mentioned  previously, 
this  will  have  the  effect  of  reducing  the  contrast  of  the  inner  details  of  these  projections, 
and  the  effect  on  the  image  is  to  eliminate  contrast  within  even  the  intersection  of  the 
correct  support  and  the  incorrect  support.  On  those  projections  that  have  support  values 
that  are  much  too  wide,  it  is  the  smoothing  terms  that  dominate.  In  order  to  lower  the 
overall  energy  of  the  sinogram  (that  is,  the  energy  term  in  the  MmIcov  random  field) ,  the 
vertical  pair-potentials  or  equivalently,  the  vertical  derivatives  should  be  small.  Therefore, 
these  projections  tend  to  become  as  smooth  as  possible  over  the  prescribed  support  and 
contribute  to  the  image  a  “shadow”  ellipse  which  corresponds  to  the  incorrect  support. 

Fig.  12  shows  a  sequence  of  reconstructions  that  have  kept  k  to  the  constant  5.0,  but 
v&iy  the  orientation  of  the  assumed  object  support.  In  these  reconstructions,  we  have  used 
the  support  of  an  ellipse  that  has  the  same  size  and  eccentricity  at  the  true  object  support, 
but  has  been  rotated  in  the  counter-clockwise  direction  by  (a)  0.0  degrees,  (b)  15.0  degrees, 
(c)  30.0  degrees,  and  (d)  45.0  degrees. 

This  set  of  experiments  shows  that  a  modest  choice  of  k,  together  with  a  less  severe 
support  error  will  produce  an  image  that  retains  many  of  the  details  of  the  true  image  with 
only  a  small  “shadow”  due  to  the  incorrect  support.  However,  it  is  clear  that  an  incorrect 
support  estimate  can  produce  results  much  worse  than  having  not  introduced  any  support 
information  whatsoever  (compare  these  results  to  that  of  Pig.  11a). 

In  Fig.  13  we  show  a  sequence  of  reconstructions  that  have  used  k  =  5.0,  but  with 
support  which  is  the  incorrect  size.  Figs.  13a  and  13b  show  two  cases  where  the  support  is 
too  small,  find  Figs.  13c  and  13d  show  two  cases  where  the  support  is  too  l2urge.  Overall, 
the  size  of  the  support  increases  from  Fig.  13a  to  Fig.  13d.  The  reconstruction  using  the 
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correct  support  and  k  =  5.0  may  be  seen  in  Fig.  12a. 

We  may  conclude  from  this  set  of  experiments  that  is  is  preferable  to  err  on  the  side 
of  using  a  support  estimate  that  is  too  large  than  too  small,  in  general.  Although  the 
boundaries  are  not  as  sharp  when  the  support  is  too  large,  the  loss  of  contrast  in  the 
interior  and  the  effect  of  double-boundaries  for  small  support  is  much  more  undesirable. 

5.5  Sparse-Angle  Studies 

Fig.  14  shows  the  results  of  several  sparse-angle  experiments.  The  (a)  and  (b)  images 
correspond  to  the  15- view  and  10- view  lOdB  cases  respectively,  where  7  =  0.05  and  ^  =  0.01 
and  the  support  is  known  and  k  =  10, 000.  The  (c)  and  (d)  images  correspond  to  the  15- 
view  and  10-view  cases,  respectively,  with  the  same  smoothing  coefficients,  but  with  k  =  0.0 
—  i.e.,  no  known  support  information  is  used. 

This  experiment  demonstrates  nicely  the  potential  of  the  algorithm.  In  either  sparse- 
angle  case,  the  contrast  of  the  image  is  improved  over  those  in  Fig.  7  dramatically.  And 
while  the  boundary  is  quite  sharp  as  expected  in  the  case  of  k  =  10, 000,  it  is  quite  clear 
what  the  shape  of  the  object  is  in  case  of  k  =  0.0.  The  loss  of  contrast  in  the  interior  of  the 
ellipse  when  k  =  10, 000  remains  evident  in  these  experiments,  however. 

5.6  Limited- Angle  Studies 

Fig.  15  shows  the  results  of  several  limited-angle  studies.  The  (a)  and  (b)  images  are 
reconstructions  obt2dned  with  known  support  (with  k  =  10, 000)  from  the  two  limited-angle 
cases.  The  experiment  resulting  in  panel  (a)  uses  the  first  40  (left-most)  projections,  whereas 
panel  (b)  uses  the  last  40  (right-most)  projections.  Panels  (c)  eind  (d)  correspond  to  the 
same  observations  as  in  (a)  and  (b),  respectively,  but  in  these  cases  no  support  information 
was  used.  As  in  the  sparse-angle  studies,  the  smoothing  coefficients  for  ail  four  studies  were 
7  =  0.05  and  jd  =  0.01. 

These  limited-angle  studies  show  behavior  which  is  similar  to  the  sparse-angle  studies. 
The  boundary  of  the  ellipse  is  quite  sharp,  as  expected,  in  the  case  of  k  =  10, 000,  and  there 
is  an  accompanying  loss  of  contrast  in  the  interior.  The  images  generated  using  k  =  0.0 
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have  different  problems,  however.  In  particular,  the  image  in  Fig.  15c  shows  good  contrast 
in  the  letters  in  the  interior  but  is  unable  to  provide  any  boundary  definition  on  the  long 
sides  of  the  ellipse.  This  is  because  the  leftmost  40  projections  which  are  observed  view  the 
ellipse  from  the  broadside,  and  as  such  do  not  contain  information  about  the  narrow  ellipse 
dimension.  The  image  in  Fig.  15d  suffers  from  the  opposite  problem.  There  is  a  loss  of 
definition  of  the  letters  in  the  interior  because  many  of  the  projections  that  would  normally 
be  obtained  from  the  broadside  of  the  ellipse  are  missing.  It  is  in  the  first  case  that  support 
knowledge  can  aid  tremendously;  unfortunately,  when  projections  from  the  broadside  of  the 
ellipse  are  missing,  there  is  little  that  our  method  can  do  to  provide  any  additional  clarity 
of  the  interior  detail. 

6  DISCUSSION 

It  is  generally  acknowledged  in  the  computed  tomography  literature  that  in  the  case  of  noisy 
and  limited-angle  or  sparse-angle  data,  prior  knowledge  is  essential  in  order  to  obtain  good 
reconstructions.  We  have  focused  on  three  types  of  prior  knowledge: 

•  Line  integrals  close  in  either  angle  or  lateral  displacement  tend  to  be  similar  in  value, 

•  Radon  tremsform  functions  zure  constrained  to  a  certain  functioned  subspace,  and 

•  Knowledge  of  the  convex  hull  of  the  object  is  equivalent  to  knowledge  of  a  particular 
region  of  support  of  the  object’s  Radon  transform. 

In  Section  3,  we  developed  a  Markov  random  field  (MRF)  model  of  sinograms  that  contains 
prior  information  about  the  meuss  and  center  of  meiss  of  the  unknown  sinogram,  the  convex 
support  of  the  object,  *md  the  expected  similarity  of  line  integrals  which  are  close  in  either 
angle  or  lateral  displacement.  The  primal-dual  MAP  estimation  algorithm  developed  in 
Section  4  is  based  on  a  variational  formulation  that  leads  to  an  efficient  solution  of  the 
original  Markov  random  field  MAP  estimation  problem.  Even  with  the  necessity  of  keep¬ 
ing  and  updating  Lagrange  multipliers,  this  method  is  fast  and  memory  efficient,  and  is 
parallelizable  in  both  the  primal  phase  and  the  dual  phase. 
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The  simulations  presented  in  Section  5  show  the  range  of  results  that  may  be  obtained 
using  this  approach.  The  improvement  over  the  unprocessed  CBP  reconstructions  is  quite 
dramatic  in  all  cases  where  the  support  was  correct  or  nearly  correct.  In  particular,  the 
boundary  of  the  ellipse  is  made  much  sharper,  and  the  letters  within  the  ellipse  can  be  made 
out  in  all  of  the  processed  cases,  and  in  none  of  the  unprocessed  cases. 

In  most  cases  the  convex  support  of  the  object  is  not  known  a  priori.  In  other  research, 
we  have  reported  several  methods  to  estimate  the  convex  hull  of  objects  from  the  available 
data  [43], [44], [32].  These  methods  require  two  steps:  1)  estimation  of  the  support  values 
from  observed  projections  and  2)  estimation  of  a  complete  set  of  feasible  support  values 
for  all  projections  (including  the  ones  corresponding  to  missing  observations).  The  feasible 
support  values  uniquely  identify  a  convex  polyhedron,  which  serves  eis  our  estimate  of  the 
convex  support  of  the  object.  In  the  second  step,  different  types  of  prior  geometric  knowledge 
about  the  shape  of  the  object  may  be  used  to  force  estimates  to  be  circular,  elliptic,  or  have 
smooth  boundaries.  For  exeunple,  knowledge  that  the  MIT  ellipse  is  an  ellipse  leads  to 
nearly  perfect  support  estimation  from  the  noisy,  sparse-angle,  and  limited-angle  cases  used 
in  Section  5,  even  though  the  size,  eccentricity,  and  orientation  of  the  ellipse  is  not  known 
a  priori.  Estimating  the  convex  support,  mass,  and  center  of  meiss  of  the  object  is  viewed 
as  a  part  of  a  hierarchical  reconstruction  algorithm  in  [32].  With  these  steps  in  place,  the 
only  user  inputs  that  eire  required  are  the  values  of  the  smoothing  coefficients  /?  and  7. 

In  this  paper  only  two  consistency  constrsdnts  were  imposed:  m&aa  and  center  of  mass. 
A  method  that  incorporates  a  much  larger  number  of  consistency  constraints  and  that 
requires  no  prior  geometric  knowledge  related  to  these  constraints  —  e.g.  mass  and  center 
of  mass  —  is  described  in  [32].  The  two  methods  contrast  in  the  following  way.  The 
approach  described  in  this  paper  characterizes,  through  the  mass  and  center  of  mass,  a 
functional  subspace  in  which  the  desired  sinogram  must  lie,  and  forces  this  to  happen  using 
Lagrange  multipliers.  The  alternate  approach  described  in  [32]  characterizes  the  subspace 
orthogonal  to  the  desired  sinogram,  and  again  uses  Lagrange  multipliers  to  achieve  this  goal. 
The  alternate  method  is  also  a  generic  primal-dual  optimization  algorithm,  however,  the 
Lagrange  multipliers  are  scalars  rather  than  functions.  The  primal  stage  is  almost  identical 
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to  that  given  herein,  but  the  dual  stage  generally  requires  more  computation  and  has  less 
potential  for  parallelism. 
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Figure  1;  The  geometry  of  the  2-D  Radon  transform. 
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Figure  4:  (a)  The  vertical  and  (b)  horizonted  boundaries  of  the  sinogram  MRF. 
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Figure  5:  The  MIT  ellipse. 
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Figure  6:  (a)  A  noise-free  sinogram,  (b)  and  its  reconstruction,  (c)  A  noisy  sinogram 
(SNR=3.0dB),  (c)  and  its  reconstruction. 
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Figure  7:  Reconstructions  from  a  noisy  sin< 
(b)  10  sparse  views,  (c)  left-most  40  views, 


n  (SNR=10.0dB)  from  (a)  15  sparse  views, 
(d)  right-most  40  views. 


Figure  8:  Estimates  produced  by  the  primal-dual  algorithm  with  (a)  7  =  0.05  and  /3  =  0.01, 
(b)  7  =  0.5  and  ^  =  0.01,  (c)  7  =  0.05  and  =  0.1,  and  (d)  7  =  0.005  and  =  0.001. 
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Figure  10:  Effect  of  known  support  for  (a)  k  =  0.0,  (b)  k  =  5.0,  (c)  k 

K  =  10,000.0. 


=  10.0,  and  (d) 
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Figure  11:  Effect  of  incorrect  support  for  (a)  k  —  0.0,  (b)  k  =  5.0,  (c)  k  —  10.0,  and  (d) 

K  =  10,000.0. 
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Figure  12:  Effect  of  using  an  incorrect  support  which  is  rotated  counter-clockwise  by  (a) 
0.0  degrees,  (b)  15,0  degrees,  (c)  30.0  degrees,  and  (d)  45.0  degrees. 
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Figure  13;  EflFect  of  using  an  incorrect  support  which  is  too  small  in  (a)  and  (b)  and  too 
large  in  (c)  and  (d). 
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Figure  14:  Spaxse-angle  studies  with  7  =  0.05  and  =  0.01.  (a)  15  observed  projections  and 
known  support,  (b)  10  observed  projections  and  known  support,  (c)  15  observed  projections 
and  no  support,  (d)  10  observed  projections  and  no  support. 
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Figure  15;  Limited-angle  studies  with  7  =  0.05  and  /?  =  0.01.  (a)  Left  40  projections  and 
known  support,  (b)  Right  40  projections  and  known  support,  (c)  Left  40  projections  and 
no  support,  (d)  Right  40  projections  and  no  support. 
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